Engineering Ising-XY spin models in a triangular lattice 
via tunable artificial gauge fields 
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Emulation of gauge fields for ultracold atoms 
provides access to a class of exotic states arising 
in strong magnetic fields. Here we report on 
the experimental realisation of tunable staggered 
gauge fields in a periodically driven triangular 
lattice. For maximal staggered magnetic fluxes, 
the doubly degenerate superfluid ground state 
breaks both a discrete Z2 (Ising) symmetry and a 
continuous U(l) symmetry. 

By measuring an Ising order parameter, we ob- 
serve a thermally driven phase transition from an 
ordered antiferromagnetic to an unordered para- 
magnetic state and textbook-like magnetisation 
curves. Both the experimental and theoretical 
analysis of the coherence properties of the ultra- 
cold gas demonstrate the strong influence of the 
Z2 symmetry onto the condensed phase. 

Phase transitions in systems with combined continuous 
and discrete symmetries are fundamentahy different from 
their purely continuous and discrete counterparts. The 
interplay between different types of excitations in the var- 
ious degrees of freedom can lead to a complex behaviour 
and coupling of the associated order parameters [1-5]. A 
paradigm example is the fully frustrated XY model on a 
triangular lattice. It combines vector spin- type symme- 
tries with discrete chiral degrees of freedom, which result 
in the famous spin-chirality coupling at low temperatures 
[6]. However, experimental studies in solid-state systems 
are challenging in view of implementing and isolating an 
XY model Hamiltonian [7-9]. 

Ultracold bosonic quantum gases in optical lattices, on 
the other hand, constitute a highly versatile system with 
an extraordinary degree of control [10, 11]. In particu- 
lar, the recent experimental realisations of artificial gauge 
potentials for bulk [12-15] and optical lattice systems [16- 
19] allow for the investigation of new physical regimes, not 
realisable in condensed matter systems. 

Here, we demonstrate the realisation of a system with 
combined U{1) and Z2 symmetries using ultracold atoms 
submitted to artificial gauge fields. Our experimental 
setup consists of an ultracold gas of ^^Rb atoms held in 



a two-dimensional triangular lattice [20] (see Fig. la). At 
each lattice site j with particle number TVj, the weakly 
interacting superfluid gas can be described by the local 
order parameter (a^) = y^V^e^"^^ . As a central aspect, 
the local phases (pj are mapped onto classical XY spins 
Sj = (cos(/pj, sin(/pj), where the tunneling matrix elements 
between neighbouring lattice sites correspond to the spin- 
spin coupling parameters. Such classical spins possess a 
continuous degree of freedom. In presence of a long-range 
order, analogous to the onset of Bose-Einstein condensa- 
tion (BEG), the order parameter assumes an arbitrary, 
but fixed phase, thus breaking the continuous U{1) sym- 
metry [21]. 

Beyond that, we experimentally engineer strong stag- 
gered gauge fields, which generate an additional discrete 
Z2 symmetry in our system. The resulting magnetic flux 
induces cyclotron-like mass currents around each plaque- 
tte. The two possible chiralities of these currents circu- 
lating around a single plaquette correspond to a discrete 
Ising-like order parameter. Furthermore, the tunability of 
the artificial gauge fields enables us to bias the Z2 order 
parameter, in analogy to a longitudinal external magnetic 
field in the Ising-spin model. 

The result is a flexible model system which allows us to 
study the temperature-dependent behaviour and interplay 
of the discrete and continuous order parameters. 

In the work presented here, the complex tunneling ma- 
trix elements, necessary to generate staggered fluxes, are 
created by accelerating the lattice potential along a closed 
orbit. A suitable periodic forcing [18, 22] results in the 
following effective Bose-Hubbard Hamiltonian: 
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where the spatial degrees of freedom perpendicular to the 
lattice have been omitted for clarity (see Supplementary 
Material). Here ^' 

operator of a boson at lattice site j, Uj - 
respective number operator, and U is an on-site repul- 
sion. In the kinetic term, the summation over the near- 
est neighbours is directional as 0ji = —6ij. The hopping 
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FIG 1. Illustration of the triangular optical lattice, the 
artificial gauge fluxes, the phase distribution, and the 
mass currents, a, The triangular optical lattice is created 
by the interference of three running- wave laser beams. In the 
experiment, roughly 2000 triangular plaquettes are occupied. 
b, Orientation of the arguments of the tunneling matrix el- 
ement around a plaquette. c, A strong artificial, staggered 
gauge field is applied to the lattice system. Crosses (dots) 
correspond to inwards (outwards) pointing gauge fluxes, d. 
An accumulated flux of ±7v around neighbouring plaquettes 
results in two, energetically degenerate, phase configurations. 
These phase configurations lead to opposite chiralities in the 
mass currents around the plaquettes. e. The current on each 
plaquette defines the orientation of an Ising-type spin. 



parameters along the directions 2^3 and 3^1 (see 
Fig. lb) are equal and denoted as \J^\e^^ in the follow- 
ing. Experimentally, \J2i\e'^^^ = \J\e'^ and \J'\e'^' can be 
tuned independently of each other. The total phase ac- 
cumulated on a closed path around one triangular pla- 
quette reflects the gauge flux through the cell, defined as 
<l> = <I>A = ((9 + 26') mod 27r = - c^^. The global accelera- 
tion of the lattice potential realised here induces fluxes 
with opposite sign for upwards and downwards pointing 
plaquettes, as depicted in Fig. Ic. 

The triangular lattice is fully frustrated for stag- 
gered fluxes of maximum magnitude tt. For this ex- 



treme case the flux structure is not unique (since 
— 7rmod27r = 7rmod27r) and the two flux patterns 
sketched in Fig. Ic are equivalent. This equivalence leads 
to two energetically degenerate vector spin configurations, 
as depicted in Fig. Id. The staggered currents induced by 
the gauge fluxes display the same degeneracy (see Fig. le). 
Note that both the argument Oij of the tunneling param- 
eter and the relative orientation {cpj —(fi) of the XY spins 
influence the mass current (jij) along one lattice bond: 
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The strong interplay between the chirality of the 
cyclotron-like mass currents (Ising parameter) and the 
XY spin long-range order induces the coupling between 
the broken Z2 and U{1) symmetries in our system. 

The presence of staggered gauge fluxes has a direct sig- 
nature in the momentum space. The single-particle dis- 
persion relation of the lattice is indeed strongly deformed: 



e{k)= -2| J|cos(k-ai -6>) 
-2|JVos(k-a2-6>') 
-2\f\cos{k-Sis-0') 
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where the a^ are the lattice directions (see Methods). For 
<l> = TT, it exhibits two degenerate minima with opposite 
ky values within the first Brillouin zone, while for fluxes of 
^ = TT ± /3 this degeneracy is lifted (see Fig. 2). For ultra- 
cold bosonic gases, the changes in the momentum space 
occupation can be easily observed with standard time-of- 
flight (TOF) imaging techniques, where the in-situ quasi- 
momentum distribution is converted into position infor- 
mation. Figure 2a shows TOF images summed over many 
experimental realisations for three amplitudes of the stag- 
gered gauge fluxes. For <!> = tt, both momentum modes 
are equally populated on average. For a strong bias flux 
of P = 0.2 TT, we externally drive the system into one of 
the minima in the flrst Brillouin zone. The correspond- 
ing dispersion relations, plotted in Fig. 2b, illustrate the 
deformations induced by the different values of the gauge 
fluxes. Figure 2d-e demonstrates the experimental control 
over the degeneracy between the two minima in the flrst 
Brillouin zone. In analogy to the effect of a longitudinal 
magnetic fleld in the Ising model, the ability of tuning the 
flux strength $ thus enables us to bias the system towards 
one of the two minima. 

The measured quasimomentum distribution contains in 
fact more information, reflecting both symmetries of the 
system. A long-range order of the XY spins, which breaks 
the U{1) symmetry, implies that the momentum distribu- 
tion is singular. 

The two possible chiralities of the mass currents cor- 
respond to quasimomenta in complementary parts of the 
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FIG 2. Effect of staggered gauge fluxes in momentum space, a, Experimentally observed occupations of the momentum 
states within the lowest Bloch band in TOF images, averaged over about 200 single-shot realisations, and b calculated dispersion 
relations for the given values of the gauge flux $. The red hexagon indicates the first Brillouin zone, c. The value of the 
real-space Ising order parameter corresponds to the occupation of a triangular mask in quasimomentum space. Atoms in the +1 
( — 1) regions correspond to positive (negative) chirality. d, A zoom into the central region along ky of the TOF images shows 
the relative occupation of the two Ising modes as a function of the gauge flux. The observation is in good agreement with the 
position and the relative importance of the minima in the band structure as shown in e. 



Brillouin zone (see Supplementary Material). Measuring 
the differential occupation in the two momentum classes, 
depicted as upwards and downwards pointing triangles in 
Fig. 2c, gives access to the mean chirality of the system. 
This analogue to the Ising magnetisation is analysed in 
the following. 

As a central result, a thermally induced phase tran- 
sition between an antiferromagnetic and a paramagnetic 
phase can be observed. Figure 3a shows a statistical anal- 
ysis of consecutive, individual experimental realisations 
for $ = TT for three different temperatures. For individ- 
ual measurements, the Ising- type magnetisation fluctu- 
ates. At the lowest temperature achieved, its statistical 
distribution clearly shows the spontaneous breaking of the 
Z2 symmetry into two individual modes. When the tem- 
perature is increased, the spontaneous magnetisation de- 
creases and finally vanishes when the system crosses the 
phase boundary to an unordered paramagnetic state. The 
simultaneous observation of both Ising states in a single 
experimental realisation is very likely due to spatial phase 
separation of different chiralities similar to the formation 
of magnetic Weiss domains. 

The bias fiux /3 impacts onto the occupation of the two 
Ising states. In Fig. 3b-c, the measurement of the mag- 
netisation as a function of the gauge fiux in the three 
temperature regimes is presented. For each value of the 



gauge fiux, the statistical distribution of the magnetisa- 
tion is represented by normalised histograms in row b. 
Row c shows the maxima of a uni- or bimodal probabil- 
ity distribution fitted to the raw data (see Supplementary 
Material). 

For a large bias fiux P the system is completely magne- 
tised in one of the two Ising states as expected for an Ising 
spin system subjected to a longitudinal magnetic field. 
Below the critical temperature and in the vicinity of fiux 
<l> = TT, we can identify two branches of favored magnetisa- 
tions which correspond to the occupation of the two Ising 
states. This behaviour cannot be explained for a system 
in thermal equilibrium. Indeed, already the presence of a 
small external magnetic field suppresses the condensation 
in the state with higher energy. However, this state cor- 
responds to a local minimum and the system can become 
metastable. The finite occupation probability of the ex- 
cited Ising state stems from non-adiabatic dynamics. The 
amplitude of the artificial gauge field is progressively in- 
creased to its final value. During this preparation ramp 
the dispersion relation becomes fiat and thus the energy 
barrier between the two states is increased from almost 
zero to the final value (see Supplementary Material for 
more details). Therefore, a finite probability exists for 
the system to be trapped in the local minimum. 

The experimentally observed metastability arises from 
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FIG 3. Measurement of the statistical distribution of the chiral magnetisation, a, The statistical distribution of the 
magnetisation (A — v)/(^ + V) obtained from consecutive single experimental realisations at flux strength $ = 7r (left) and the 
corresponding histograms (right) are shown for three different temperatures, b, Histograms revealing the statistical distribution 
of the sample magnetisation for different temperatures and gauge fluxes. For each of these histograms about 200 individual 
measurements have been recorded. The colour code corresponds to the normalised amplitudes of the histograms, c, Maxima of 
Gaussian probability distributions which are fitted to the raw data. For bimodal distributions, the point size represents their 
relative weighting. 



the repulsive interactions between the atoms, which pre- 
vents fragmentation of the state. This is supported by the- 
oretical calculations including these interactions. Namely, 
the free energy of the system can be evaluated up to the 
first-order correction in the interaction strength [23] (see 
Supplementary Material). At low temperatures, the ef- 
fective free energy develops two minima. Condensation in 
one of the two minima is equivalent to the spontaneous 
breaking of the Z2 symmetry. The energy scale protecting 
the metastable minimum is the mean-field energy, which 
is large compared to the temperature. This conclusion is 



also confirmed by a Bogoliubov-de Gennes theory, includ- 
ing the second-order correction with respect to interac- 
tions (i.e. the quantum fluctuations) (see Methods). For 
this reason, we can observe the metastable states. 

At higher temperatures the entropic contribution to 
the free energy merges the two minima into one. There- 
fore, no metastable state is expected above the Z2 criti- 
cal point. The measured occupation of metastable states 
with a magnetisation opposite to the bias field is therefore 
a non-equilibrium signature of the phase transition. 

The equilibrium state of the studied three-dimensional 
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FIG 4. Magnetisation curves obtained via a classical 
Monte Carlo simulation. At the lowest temperature (T), a 
non-zero spontaneous magnetisation is well reproduced, which 
disappears as the temperature is increased {T' < T"). Tising 
denotes the critical temperature for the Z2 symmetry breaking. 
The apparent scattering reflects the fluctuations of a series of 
MC simulations. 



system can be investigated via a classical Monte Carlo 
(MC) approach. Here, all relevant experimental parame- 
ters, including the overall confinement, have been taken 
into account in the simulation (see Supplementary Mate- 
rial). Figure 4 shows the magnetisation curves for three 
different temperatures, generated by extracting the chi- 
rality from the momentum distribution. For flux <l> = tt, 
the thermally driven phase transition from an ordered 
state showing spontaneous magnetisation to an unordered 
state is reproduced, and overall similar to the experimen- 
tal data. No finite occupation of the metastable minimum 
can be observed, since the MC simulation generates the 
ground state of the system. 

While the statistical distribution of the magnetisation 
quantifies the Z2 symmetry breaking, the sharpness of the 
momentum peaks is a measure for the long-range phase 
coherence connected to the U{1) symmetry. The peak full- 
width-half-maximum (FWHM), extracted from the exper- 
imental data presented in Fig. 3b-c, is shown in Figure 
5a. As expected, for each flux value the peak width in- 
creases with the temperature, monitoring the decreasing 
long-range order. More remarkably the U{1) order param- 
eter depends strongly on the gauge flux strength. For a 
deeper understanding of this behaviour, the critical tem- 
perature for BEC {Tc) has been theoretically evaluated in 
the weak-coupling approximation of the free energy and is 
plotted in Fig. 5b (see Supplementary Material). For mea- 
surements realised at a fixed temperature, the coherence 
length of the gas should decrease in the vicinity of <l> = tt, 
where Tc displays a pronounced cusp. The measured full- 




FIG 5. Experimental and theoretical evaluations, re- 
lated to the U{1) order parameter, a, Measured FWHM 
of the momentum peaks for the three different temperatures 
Ti < T2 < T3 . The width of the momentum peaks is a mea- 
sure for the loss of long-range coherence of the sample, b, Free 
energy in units of J^^ / iim^ as a function of flux and temper- 
ature. The condensation transition (red line) shows a cusp at 
TT-flux. 



width-half-maximum is limited by the finite time-of- flight, 
but the short coherence lengths expected in the vicinity of 
$ = TT are nicely reproduced. The observed increase of the 
measured FWHM symmetric to <l> = tt is in good agree- 
ment with the theory. Similar conclusions follow from the 
exact thermodynamic analysis of the non-interacting gas 
(see Supplementary Material). 

In conclusion, we have realised a model system with 
Ising-type Z2 and global /7(1) phase symmetry by apply- 
ing strong gauge fields to bosonic atoms in a triangular 
optical lattice. For classical two-dimensional XY systems 
with coupled spin and chirality degrees of freedom, theory 
predicts that the system first breaks the Z2 chiral sym- 
metry and then the U(l) symmetry as the temperature 
is reduced [24]. However, the exact nature of these phase 
transitions, which are strongly linked by combined excita- 



tions, has long been debated [25, 26]. Only recently, pre- 
cise Monte Carlo simulations could resolve the two tran- 
sitions and identify their universality classes [27, 28]. The 
analysis of the coherence properties of the 3D ultracold 
gas demonstrates the strong influence of the Z2 symme- 
try breaking onto the BEC phase, revealed as a drastic 
reduction of the coherence length. In future, it will be in- 
teresting to investigate the coupling between these phase 
transitions and its influence on their critical behaviour, 
which is however experimentally challenging. In addition, 
the occupation of metastable states with a magnetisation 
opposite to the bias field is a non-equilibrium signature of 
the Ising-like phase transition. This constitutes a funda- 
mental, defining property of such phase transitions, which 
is observed here in the field of ultracold atoms. 

This work paves the way to further studies of artificial 
magnetic properties of ultracold quantum gases in optical 
lattices. Combinations of the two-dimensional control of 
the complex tunneling parameters reported here with su- 
perlattices [29] in different lattice geometries (triangular, 
hexagonal, or kagome) promise to give deeper insights 
into a variety of magnetic systems [5, 11]. 
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Computing (NIC) for providing us with computing time 
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METHODS 

The triangular optical lattice. The two dimensional, 
triangular optical lattice as depicted in Fig. Id is created by three 
running laser beams with actively stabilised phases that intersect 
in the xy-plane at an angle of 120°. The beams are derived from 
a Ti:sapphire laser at wavelength Ai^=830nm, creating a 2D 
lattice potential y(r) = — Vb Xli <^<^s(t>2i') with a lattice spacing 
of (i = 2A/,/3 = 533nm. The reciprocal lattice directions are 
bi= 6/2(1, v^,0), b2 = 6(l,0,0) and ba = V2(-l, v^, 0), where 
6 = 27rv^/Ai,, corresponding to the real-space lattice directions 
ai = d(0, 1, 0), a2 = (i/2(v^, -1, 0) and as = -d/2(v^, 1, 0). 

Experimental preparation. We create Bose-Einstein conden- 
sates of (1.5 — 2.5) X 10^ ^^Rb atoms in a crossed optical dipole trap. 
Within 100 ms, we subsequently ramp up the optical lattice to a final 
lattice depth of (4.6 ±0.1) Erec (Erec = hx 3.33 kHz) which leads to 



mately 2100 to 2600 elongated tubes with a mean occupancy in the 
range of 70 to 95 atoms (175 to 235 in the center). The temperature 
of the system is increased by holding the atoms longer in the lattice 
before applying the artificial gauge fluxes. 

Note that the system under study is three-dimensional. There- 
fore we observe a BEC transition instead of a Kosterlitz Thouless 
transition as expected for a pure two-dimensional system. 

Lattice shaking. Staggered fluxes in the triangular lattice are 
realised by a global periodic motion of the optical lattice around a 
closed orbit x(t) = — Ax cos{Cjt)ex — Ay [sm{u;t) -\- 6 sin(2ujt) / 4] e^, 
with Co ^ 27T X 2.8 kHz. The amplitude of the staggered flux can be 
accessed by the control parameter 6. In the reference frame of the 
moving lattice, this results in a force 

F(t) = —m5t(t) = —Fx cos(Cut)ex — Fy [s\n{Cjt) + 5s\n{2Cbt)] e^ (5) 

acting on the atoms. Experimentally, the trajectory is re- 
alised by modulating two of the three lattice laser beams 
with ^2/3 = ± ^cc sin(a;t) + Uy [cos(a;t) + 5 cos(2a;t)/2] , where 

Ux = Co Ax / {V^d) ^ Vy =CbAy/d and Ax,y = Fx ,y / {mCo'^) . Hereby, the 
shaking amplitudes Ux and Uy are linearly increased to their flnal 
values in 50 ms after the condensate is loaded into the initially 
resting lattice. Time-averaging the projection of the force onto the 
bonds of the elementary plaquette now leads to a renormalisation 
of the tunneling matrix elements J^^^^ — )► J^^ in the xy-plane [30]. 
The absolute values of the effective tunneling matrix elements are 
on the order of \J^^\ ^ 0.4 J^^^^. 

Symmetries of the system. With U / J^^ = 1.2 — 1.4 and 
large fllling factors, the system remains in the weakly interact- 
ing regime. Here, the local wavefunction has well deflned phases 
(^j on each lattice site, which correspond to the XY vector spins 
Sj = (cos (y9 j , sin (y9 j ) . In the case where the arguments of the effec- 
tive hopping are equal to tt (i.e. 5 = 0), the total kinetic energy of 
the atomic ensemble along the lattice directions can be written as 

= ^|J^|s,.s,. (6) 

(id) 

Equation (6) is invariant under both a discrete Z2 transformation 
and a global U(l) rotation: 
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[cos{(pj), — sm{ipj)) Z2 

( cos{(pj + ly), sin((^j + ly)) U{1) 



(7) 
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:4 X 10- 



a bare single-particle tunneling parameter of J 
As the system is only weakly conflned in z-direction with respect 
to the lattice potential - the overall external harmonic conflnement 
is WTot =277 X (31, 53, 40) Hz - the atoms form an array of approxi- 



On the contrary, the chirality changes its sign under the discrete 
transformation. The summation of this quantity over the lattice 
plaquettes corresponds to the magnetisation of the system. 

Detection and data analysis. All the information about 
the momentum distribution of the superfluid are retrieved 
from absorption images taken after 32 ms time-of-flight. The 
chirality is deflned by the spins at the corners of one elemen- 
tary plaquette as x = sgn[s2 x si + S3 x S2 + si x S3], where 
s^ X Sj = Si^xSj,y — Si^ySj^x- It can be converted to a mask in 
quasimomentum space (see Fig. 2c). By weighting each absorption 
image with this mask we obtain the total magnetisation of the 
system as shown in Fig. 3b. 

Bogoliubov theory for the metastable condensate. In a 

translational invariant system with a Bose-condensate in one of the 
two local minima of the free dispersion relation ^(k), one can add 



quantum and thermal fluctuations within Bogohubov theory. One 
obtains the quasiparticle dispersion relation 



cj(q) : 



^(q) - ^(-q) 



6(q)+6(-q) 
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for momenta q = k — kg relative to the condensate momentum kg, 
where ^(q) = ^(ko + q) — ^(ko) and with g-^^ and p-^^ denoting 
the interaction parameter and the density in the tubes respectively. 
A thermodynamic instability is indicated when a;(q) assumes neg- 
ative values for some q. For vanishing interaction Pi^Pid = 0, 
one has oj(q) = e(q) and the system is thermodynamically unstable 
as soon as the condensate is not prepared in the global minimum 
of the dispersion relation. However, flnite interactions gn^Pn^ > 
can stabilise a condensate in the upper local minimum of the disper- 
sion with a;(q) > 0. The very same mechanism leads to spontaneous 
symmetry breaking for <l> = tt by disfavoring a fractionalised conden- 
sation. The existence of a metastable state is thus directly linked to 
spontaneous symmetry breaking. 
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Supplementary material 



LATTICE SHAKING 

As stated in [SI], staggered fluxes in triangular lattices 
can be realised by a global periodic motion of the optical 
lattice around a closed orbit. The trajectory used for the 
experiments presented in this article is given by 



x(t) = — Ax cos{ujt)ex 

— Ay [sin{ujt) + ^sin(2a;t)/4] e^, 



(SI) 



where a; = 27r/T = 27r x 2.791kHz with x(t)=x(t + T). 
The important control parameter for the staggered flux 
strength is S. For ^ = all tunneling matrix elements are 
real valued and only flux strengths which are zero or tt 
can be achieved. The inertial force acting on the atoms 
in the reference frame of the moving lattice is 

F{t) = — Fx cos{ujt)ex — Fy [sm{ujt) + Ssm{2ujt)] e^, 

(S2) 
where the connection to the trajectory (eqn. SI) is given 
by Ax = Fx / {muj'^) and Ay = Fy / {muj'^) . Experimentally 
the forcing of the atoms in the lattice is realised by fre- 
quency modulating two of the three lattice laser beams 
with 



Ai^2/3 = ±^x sin(a;t) 

-{-Uy [cos{ujt) + S cos{2ujt)/2] , 



(S3) 



where Ux = Fx / {V^d'f^^) and Uy = Fy / {dmuj) . The 
renormalised tunneling matrix elements due to the time 
averaging over one cycle are 

Tbare pT 

-j^ j dtexp(iW-i,(i)/^), (S4) 



7eff . 




FIG S2. Minima in the first Brillouin zone in presence 
of gauge fluxes, k-space separation between the ground- and 
excited state quasi momenta in dependence of the anisotropy 
parameter a and the staggered flux strength $. 



with 



Wi 



i{t)=-f dTF{T)aj. 

J —OO 



(S5) 



The vectors slj describe a closed path around one elemen- 
tary plaquette of the lattice. Taking advantage of the 
symmetries in our system, the tunneling matrix elements 
are written as J = J21 and J' = J32 = J13 in the following. 
Fig. Sla depicts the numerical solutions for the Peierls 
phases and the resulting staggered flux according to the 
equation (S4). In Fig. Sib the magnitude of the different 
effective tunneling matrix elements are shown. Note that 
the difference between the magnitude of the tunneling ma- 
trix elements is on the order of a few percents. Therefore, 
it only has a weak influence on the dispersion, as will be 
detailed in the next section. 
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FIG SI. Tunability of the complex tunneling parame- 
ters, a, The Peierls phases along the bonds and the resulting 
flux strength through an elementary plaquette are plotted as 
a function of the control parameter S. b, Magnitudes of the 
three effective tunneUng matrix elements in units of the bare 
tunneling amplitude. 



THE DISPERSION RELATION 

The lowest band dispersion relation for the triangular 
lattice in presence of complex tunneling matrix elements 
is described by 



s{k)= -2\J\cos{dky-02i) 

- 2\f\ cos {d [Vskx - ky] /2 - 632) 

-2\f\cos{d[V3kx^ky]/2- 



(S6) 



ns 



The gauge invariant quantity of the system is the stag- 
gered flux strength $. On the contrary, the hopping ar- 
guments Oij depend on the chosen gauge. A change of 
gauge corresponds to a translation of the band structure. 
The specific choice of the gauge ^21 = tt + /3, 632 = tt and 



^13 = TT yields the simplified expression for the dispersion 
relation 



e{k) = ^2\ J \ cos {dky-/3) 

+4| J'l cos {dk^V3/2) cos {dky/2) . 



(S7) 



The x-component of the quasi-momentum for the local 
minimum is qmin,x = 2ii /{d\/i). For <l> = tt the correspond- 
ing ^-component can be written as: 



Qmin^y 



_ for Of > 2 

I ± I arccos ( ^ ) for a < 2 



(S8) 



where a — J j J' is the anisotropy parameter of the lattice. 
The numerical results for the quasi-momentum separation 
between ground-state minimum and metastable minimum 
are shown in Fig. S2. Since in our case the anisotropy 
parameter remains close to unity, the quasi-momentum 
separation is only weakly depending on the flux. 



SWITCHING THE GAUGE FLUXES 

After the lattice potential has been ramped to its fi- 
nal depth of 4.6£^rec7 the tunneling matrix elements are 
all real and positive valued. In order to induce non-zero 
gauge fluxes, the frequency modulation of the laser beams 
is slowly turned on by increasing the frequency amplitudes 
Vx and Vy linearly over a time Tr = 50 ms. Depending 
on the value of the control parameter ^, staggered gauge 
fluxes with different final amplitude can be realised. 

Since the ramping time scale is slow compared to the 
orbital motion of frequency a; = 27r x 2.791kHz, the sys- 
tem has well defined tunneling matrix elements during the 
switching procedure. The time resolved evolution of the 
phases and amplitudes of the effective hopping elements 
during the ramping of v^ and Vy are shown in Fig. S3. It is 
important to note that the fluxes are rapidly switched to 
their final amplitude. This corresponds to a quench into 
the final state and explains the non-adiabatic behaviour 
described in the main text. As observed experimentally, 
slow ramps reduce the excitations in the system but since 
the absolute tunneling values become small during the 
ramp (see Fig. S3b and c), the process cannot be fully 
adiabatic in experimentally accessible time scales. 

The initial temperature of the system has been varied 
by holding the atomic sample in the lattice prior to intro- 
ducing the staggered gauge fluxes by shaking. The chosen 
durations were ms, 80 ms and 160 ms respectively for the 
three investigated regimes. 



TUBE PARAMETERS 

In order to understand the physical properties of the 
complete system in all three dimensions, it is important 




FIG S3. Switching of the gauge fluxes, a, Amplitudes of 
the gauge flux strength and b, c of the tunneling matrix ele- 
ments along the bonds 1^2 and 2 ^ 3 (3 ^ 1) respectively 
are plotted for different values of the control parameter ^ as a 
function of the time during the linear ramping of Vx and Vy. 
The time is expressed in units of the ramp time Tr. The in- 
set in c is a zoom into the region of small absolute tunneling 
matrix elements J32/i3- 



to derive some basic parameters for the array of elongated 
tubes in z-direction that are formed by the presence of the 
optical lattice in the x?/-plane. The basic ansatz for the 
wavefunction of the system is 



V(r): 



\ciWi(x,y)Qi[z), 



(S9) 



where Wi{x^y) is the single particle Wannier function of 
the 2D lattice, Ci{z) is an interaction broadened function 
along the tubes and Ni = \ci\'^ is the number of atoms in 
the tube of lattice site i. Neglecting the kinetic energy of 
the system, this leads to a Gross-Pitaevskii equation for 
the functions Ci(^)- The squared modulus of the function 
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is given by 



^i-m (ujIRI^ + ujlRly + uj^.z^) /2 



(SIO) 



where the ^x,y,z denote the overah external harmonic 
confinement and Ri^x-,Ri,y are the x^y components of 
the lattice vector of site i. g = g j dxdy\wi{x^y)\^ is the 
renormalised interaction parameter with the bare three 
dimensional interaction parameter g = AT:h^as/m. For 
^^Rb: aF = 2= + (100.4 ± 0.1)ao [S2]. A numerical cal- 
culation of the 2D-Wannier functions yields the result 
g/g = 17.2 x /im~^ for 4.6Erec- The length of the tubes 
{2ztf) is determined by the Thomas- Fermi boundaries in 
2;-direction 



ZTF-- 



l2^-m{ulRl^ 



-ujlRl. 



/2 



(Sll) 



and the number of particles in a single tube 
is Ni = 2muj1z'^Y/{'^9)' ^Y making the contin- 

uum approximation ^^ Ni -^ ^uc / ^xdyN{x^ y) with 
Ri,xi Ri,y -^ ^5 y^ the chemical potential can be calculated 
analytically: 



/^ = 



flbgAuc 3/2 ,. 



2/5 



(S12) 



where Auc = \/3/4 (P is the area of the unit cell and 
A^Tot the total particle number of the system. The energy 
scale associated with each tube i can be calculated as 
Ui=g J dz\(i{z)\^ = 3g/{bzTF)- Relevant system parame- 
ters for total particle numbers of 1.5 x 10^ and 2.5 x 10^ 
are depicted in Tab. SI. 



iVTot 




1.5 x 10^ 


2.5 X 10^ 


Asites 




2157 


2629 


Axot/Asites 




70 


95 


Anaax 




174 


237 


U/J= E^(f^^^0/(ATot0.4 J^^^^) 


1.4 


1.2 


^Tube = X)i(2 2;i,TF 


AO/Axot 


22.7/im 


25.2 /im 


PlI^=E^[N!/{2z 


'i,TF)] /Axot 


4.6 /xm"^ 


5.6 /im~"^ 



TAB SI. Calculated system parameters for different to- 
tal particle numbers. Parameters are: the number of occu- 
pied sites Asites, mean tube occupancy ATot/Asites, maximum 
tube occupancy Amax, occupation weighted ratio U/J, occu- 
pation weighted tube length /xube and occupation weighted 
ID-Density piD- 



where wq is the Fourier-transformed Wannier function 
driven by the shaking. The expectation value in the sum 
describes the coherence properties of the sample. In order 
to keep the same Wannier envelope position, the switch- 
off time is chosen to occur at the same time within one 
period for all measurements. The center is positioned in 
between the two degenerate minima (for ^ = tt) of the dis- 
persion relation at kc = {-\-27T/\^d^ 0). However, a small 
displacement remains towards ky <0. For the given en- 
velope size the displacement leads to a slightly favored 
weighting of the negative magnetisation in the performed 
measurements. This is the reason for the negative offset 
of the data presented in Fig. 3. 



TIME-OF-FLIGHT MEASUREMENTS 

After rapidly switching off all trapping potentials and 
letting the atoms fall in free space for 32 ms we take 
an absorption image of the cloud with a magnification 
of approximately 3. With this standard time-of-flight 
method, the quasi- momentum distribution of the atoms 
in the lattice can be revealed. In Fig. S4, samples of av- 
eraged time-of-flight images are shown in dependence of 
flux strength and temperature in the lattice. As stated 
before, the physics of the shaken system is described by 
a time-averaged effective Hamiltonian with renormalised 
tunneling matrix elements. The quasi- momentum distri- 
bution describing the effective model is static, while the 
only effect of the fast periodic acceleration is an overall 
oscillating envelope on top of this quasi-momentum dis- 
tribution. The density distribution in the far-field regime 
after time-of-flight is given by 



n(k): 



^0 k 



mr 



^e^k(R.- 



R, 



aja. 



(S13) 



GAUGE-INDEPENDENT CHIRALITY MASKS 

In order to determine the magnetisation of the system 
from the time-of-flight images, we introduce the chirality 
of the system 



X = sgn [S2 X Si + S3 X S2 + Si X S3] , 



(S14) 



where the scalar cross product of our two dimensional 
vector spins is defined as 



07 X. So 



e^-'s.s,- 



= Si 



Si.ySjx- 



(S15) 



The spins Si, S2, S3 are arranged clockwise around a trian- 
gular plaquette. For staggered mass currents this results 
in the same value of x foi" the two types of triangular 
plaquettes (upwards- and downwards pointing triangles 
in x-direction). The chirality can be converted into the 
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FIG S4. Quasi-momentum distributions, a, Dispersion relation £(k) for selected values of the flux strength $. b, The 
corresponding averaged time-of- flight images for the three different initial temperatures Ti, T2 and T3 show the characteristic 
population of the minima in the dispersion relation. The dashed lines in the first row of images are a guide to the eye. 
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FIG S5. Gauge dependency of the chirality masks. The 

chirality masks xg obtained from the sign of equation (S18) 
are plotted for different gauges: a ^ = O.Ott, 0' = 0.957r and b 
6 = I.Itt, 6' — I.OStt. The first Brillouin zone is indicated by a 
red hexagon. 



reciprocal space, resulting in a mask for the TOF images 

3 



X(k)=sgn 



y sin(k • 



a, 



=1 



(SI6) 



Weighting the TOF images of n(k) with this mask 
gives access to the mean chirality or, in the Ising 
picture, the staggered magnetisation of the system 
M= Jx(k)7i(k)d^k. In principle this quantity is not ex- 
act since it is gauge dependent. An observable which char- 



acterises the Z2 order parameter in a gauge independent 
way is given by the total staggered flux jTxot, that is 



A^si 



(jTot)= Yl Eo-.) 



j=i 



M 

h 



(S17) 



^(n(k))A'(k,^,^0. 



Here, n (k) is the density in momentum space and (jia ) 
denotes the expectation value of the mass current from 
lattice site j to the nearest neighbour in the direction a^. 
With the different and 0' , the specific gauge that was 
chosen is described by the weighting function 

A'(k,6>,6>0 = sin(k-ai -6>) 

+ sin(k-a2-6>0 (S18) 

+ sin(k-a3 - 0') . 

It can be used to define a set of gauge-independent chiral- 
ity masks XG(k, 6>,6>') =sgn [A' (k, 6>, 6>')]. The resulting 
masks for the specific gauge of ^ = 0.97r, 0' = 0.957r and 
^ = I.Itt, ^' = 1.057r, that correspond to flux strengths of 
$ = O.Stt and $ = 1.27r, respectively, are shown in Fig. S5. 
For simplicity, we use the mask defined by equation (S16), 
where = 0' = t:^ and Xg = Xi since the difference of the 
magnetisation data generated with gauge independent 
masks turns out to be negligibly small. 



12 




100 200 
Measurement 



25 
Occurrence 



FIG S6. Analysis of the statistical distribution of the 
magnetisation. Illustration of the bimodal fluctuation of the 
magnetisation due to the breaking of the Z2 symmetry for the 
case of low temperatures and flux $ = tt as seen in Fig. 3a. The 
solid lines on the right represent the fitted uni- (/cm = 1) and 
bimodal (/cm = 2) Gaussian probability distributions. 



STATISTICAL DATA ANALYSIS 



width and amplitude, thus hinting at a remaining dis- 
crepancy in the evaluation of the respective fits). This 
is a decisive evidence of the phase transition from an or- 
dered, ferromagnetic to an unordered, paramagnetic state. 
The resulting Gaussian density distributions (Fig. S7c) are 
in good agreement with the statistical representation of 
the data in Fig. 3b. In Fig. 3c (and again in the inset of 
Fig. S7c) the maxima of the obtained Gaussian distribu- 
tions are plotted. In the case of bimodal distributions, the 
point size represents the ratio of amplitudes of the respec- 
tive Gaussian, emphasising the smaller population of the 
metastable minimum for the measurement at Ti. 

Another indication for the disappearance of sponta- 
neous symmetry breaking for larger temperatures is the 
behaviour of the variance of the magnetisation measure- 
ments as shown in Fig. S8. Here, a notable rise of fluctua- 
tions for fluxes close to tt is evident for Ti, and, although 
less distinct, for T2, while the fluctuations remain con- 
stantly small for T3. 



As described in the main text, the spontaneous breaking 
of the Z2 symmetry manifests itself in characteristic shot- 
to-shot fluctuations of the measured total magnetisation 
of the system. Therefore, a statistical analysis of the data 
is essential in order to extract reliable information about 
properties of the raw data distribution plotted in Fig. S7a. 

For this purpose, we fit a one-dimensional Gaussian 
probability distribution with /cm = 1,2 modes to the mag- 
netisation data for each flux value. With such a soft clus- 
tering method, the actual number and properties of modes 
in the ID-distributions can be determined by comparing 
the Schwarz-Bayes criterion (SBC) for the respective fits 
[S3]. Fig. S6 illustrates a uni- and bimodal Gaussian dis- 
tribution for the case of flux <l> = tt for the measurement 
with lowest temperature, where SBCi > SBC2 favors the 
bimodal model. In order to assure reliable results, each 
fit is replicated ten times with random starting parame- 
ters, selecting the most likely output. Furthermore, the 
obtained parameters are averaged ten times so that de- 
viations due to the randomness of the initial fitting pa- 
rameters can be ruled out. In Fig. S7b the differences 
of the SBC for uni- and bimodal fits are plotted. For 
cases where SBCi < SBC2 (SBCi > SBC2) a unimodal 
(bimodal) model is favored. Note that multimodal dis- 
tributions with kM > 3 are not considered here since they 
have proven to be always less favorable as compared to 
the cases /cm = 1 and /cm = 2. 

Spontaneous symmetry breaking for fluxes close to tt is 
clearly indicated for the temperatures Ti and T2 by favor- 
ing bimodal probability distributions. On the contrary, 
no symmetry breaking can be observed for the tempera- 
ture T3 as the unimodal fit is always favored (in spite of 
the outlier for <l> = 0.9 tt, where the SBC are nearly equal 
and the two resulting Gaussian modes strongly differ in 



FREE ENERGY 

In this section we discuss the thermodynamic behaviour 
of the frustrated lattice system. We use a weak-coupling 
approximation of the free energy, which contains the non- 
interacting contribution and the first-order term, as dis- 
cussed in Ref. [S4]. 



Symmetric case 

We consider the three-dimensional lattice dispersion in 
the fully frustrated case, for $ = tt and | J| = | J'|: 



s{k) = ^2\ J\ cos {dky) 

-^2\ J\ cos {d[V3k^ 
+2|J|cos((i[V3fc^ 



ky]/2) 

ky]/2)^kl/{2m). 



(S19) 



In order to map the system on a weakly interacting 3D 
Bose system, we expand the dispersion around the two 
minima ko,A/B to second order, and write the curvature as 
an effective mass 



e{^i) 



£2 £2 c2 

2mx 2my 2mz 



(S20) 



where ^i = k — ko,i, and i = A, B with A and B denoting 
the two distinct minima in the dispersion. The in-plane 
masses are mx=my=mj = 2Ti? /{?>d?\J\)^ while the mass 
along the tube is simply the bare mass, mz=m. The 
effective 3D density is related to the ID density in the 
tubes by 77,3^3 =n^^2/{\^(P). In analogy to the isotropic 
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FIG S7. RaAv data, information criteria and fit results of the statistical distributions, a, Single shot measurements of 
the magnetisation as a function of the flux strength for the three different initial temperatures, b, Differences of Schwarz-Bayes 
criteria for a uni- and bimodal Gaussian probability distribution are plotted for each given flux value. As indicated by empty 
circles, for some cases the bimodal fit fails to converge and a unimodal fit has to be assumed as the best model, c, Resulting 
probability fits showing good agreement with the histograms from Fig. 3b. The extracted maxima of the probability distribution 
(see Fig. 3c) are shown in the insets. 
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FIG S8. Variances. Fluctuation of the magnetisation for 
the three temperatures demonstrating the symmetry breaking 
for fluxes close to $ = 7r. Each data point corresponds to the 
variance of the ID-magnetisation data for the respective flux 
value. The solid lines are Gaussian fits to the data. 



3D case, we define the two thermal wavelengths 



A^ = h/^2T:mkBT (S21) 

\j = h/y^27TmjkBT (S22) 

and A = (A^Aj)^/^. We now consider a thermal distribu- 
tion of non-interacting bosons. In analogy to the regular 
Bose gas we find for the density of excited states in the 
minima A and B 

^e,i = -3^3/2(21) (S23) 

A 

with the Bose function ^^(z) = Yll^i^^ /^^ ^^^ ^^^ cor- 
responding fugacities z±= ex.p{/jj±/kBT). The chemi- 
cal potentials /ii control the densities n^ in each mini- 
mum. The free energy of the non-interacting system is 

Aq = Ao,A + ^o,B, where 



V 



knT 



ff 95/2(^1) + UikBTlnZi if Zi < 1 



ksT 



^5/2(1) 



if Zi = 1 
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FIG S9. Free Energy in dependence of the density im- 
balance, a, Free energy A/V — {Aq + Ai)/V per volume for 
the parameters of the experiment, given in the text, as a func- 
tion of the density Ua of particles near minimum A of the 
dispersion for the temperatures Ti = 65 J^^, T2 = 67 J^^ and 
T3 = 69 J^^. We keep the total density rir^^ —tia + tib fixed and 
the energy A{nA — nr^^/2^T) is set to zero, b. Contour plot of 
the free energy A/V (in units of J^^//im^) for a wider tem- 
perature range. The two new minima appear symmetrically 
around nA — Ur^^l^^ as the temperature is lowered. 



In order to account for the interaction, we include the first- 
order term in the effective 3D interaction strength g^^^ 
which is related to the ID interaction g^^ in the tubes by 
^3D =^id2<^^/a/3- As discussed in Ref. [S4], the first order 
correction to AjV is 

^ = ^ [2(n^ + n^f - nl^ - nl^\ (S24) 

where no,A and no,B are the condensate densities in 
minimum A and B, respectively. In Fig. S9, we plot 
the free energy per volume A/F = (Aq + Ai)/y, for 
J = J®^ = /cb X 0.26nK, and for a fixed total density 
of n3j3 = 17/im~^, corresponding to a ID density of 
n^j3=6/im~^. With a ID interaction strength of 
^iD — 23.4 J^^/j.m, this results in an effective 3D interac- 
tion strength of g^^^ =6AJ^^/im^. As the temperature T 



is lowered, the free energy develops two minima symmet- 
rically around rii = n^^/2^ indicating the onset of sponta- 
neous breaking of a Z2 symmetry. Furthermore, we see 
that within this approximation the free energy barrier is 
of the order of g^^rioNo^ where no and Nq denote the den- 
sity and the number of condensed particles, respectively. 
When the condensate fraction approaches 1, the energy 
barrier per particle becomes ^sd^o^^sd^sd which is of 
the order of ks x 28 nK or 108 J®^. Since this energy is 
large compared to the temperature estimates of the exper- 
iment, it can protect the metastable states that are seen 
following the quench. We also note that in this estimate 
the breaking of the Z2 and the U{1) symmetry occur at 
the same temperature, because it is the condensate frac- 
tion that is responsible for generating two minima in the 
free energy. 



Biased case 

We now consider the case where the minima of the 
dispersion relation are not degenerate, but have an en- 
ergy difference of A = £:(ko,B) — ^(ko,A) resulting from 
a fiux value different from <l> = tt. An approximate re- 
lation between tilt energy A and flux strength $ is 
A = 10.5 J®^ X (^/tt — 1). We choose the energy minima of 
the dispersion such that 6:(ko,A) = 1^1 and 6:(ko,B) = for 
A < 0, and £:(ko,A) = and 6:(ko,B) = A for A > 0. Using 
the same approximation as in the previous section we find 
the following expression for the free energy 



A/V 



^[^5/2(Za)+^5/2(Zb)] 

^ksT [ua In(ZA) + Ub In(ZB)] 



(S25) 



+ |[2(nA + nB)' 



^0,A ■ 



^0,E 



If only the density 77,3^3 = tia + ^b is given, as it is the case 
for the experiment, only one of the density fractions, no, a 
or no,B, can be non-zero. As is apparent from equation 
(S25), the system can always lower its energy by condens- 
ing all atoms into only one of the two minima. 

We first hold the individual densities fixed to demon- 
strate the behaviour of the free energy described in the 
main text. In Fig. SIO we show the free energy A/V as a 
function of nA, with nA + nB held fixed, for different fiux 
strengths ^. In Fig. SlOa, a local minimum persists for a 
finite tilt. Here the temperature is low enough, that the 
system is condensed at finite tilt of the system away from 
^ = 7T. In Fig. SlOb we choose a higher temperature, re- 
sulting in only one global minimum being present for any 
tilt. We now only hold the total density nA+nB fixed so the 
fugacities are given by Zi = z and Z2 = zexp( — |A|//cbT), 
for A > 0, and by Zi = zexp(-|A|//cBT) and Z2 = z for 
A < 0. The density of excited states is related to the 
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FIG SIO. Free Energy behaviour in dependence of the 
flux. Free energy A/V per volume as a function of the den- 
sity Ua for various flux strengths $ and a temperature of 
a, T = 66 J^^ and b, T = 70 J^^. We keep the total density 



--nA-\-nB fixed. The energy A{nA- 



,/2,A) is set to 



zero. For the smaller temperature in a, a local minimum per- 
sists for finite tilt energy, indicating the two degenerate minima 
that exist for the symmetric case. For the higher temperature 
in b only one minimum can be seen, indicating that the system 
is supercritical. 



fugacity z through 



ne.A + n, 



e,A T /^e,B 



53/2(z)+53/2(ze-l^l/'=-^)l/A^ (S26) 



The free energy is then given by 



A/V 



keT 



55/2(2) +55/2(2 e' 



-^kBTijij^ + tib) In z — ni A 



-|A|//cbTA 



(S27) 



where i = B if A > 0, and i = A if A < 0. In Fig. Sll 
we show the free energy for the same parameters as in 
the previous section, in dependence of the temperature T 
and the tilt energy A, while keeping only the total density 
77,3^3 fixed rather than the individual densities. Note that 
here we span a much larger parameter range of the tilt 
energy than is experimentally accessible (compare Fig. 4 
in the main text). For each temperature, A(A = 0, T) has 
been set to zero. If the system is above the critical point, 
the free energy increases as the tilt energy is varied away 
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FIG Sll. Free energy A/V as a function of temperature 
and tilt energy. Here, only the total density n^^ = tia + ns 
is held constant. Above the critical point, tilting the disper- 
sion increases both the free energy (plotted here in units of 
J^^//im^), and the phase space density in the lower minimum, 
leading to condensation. While the phase transition is depicted 
by the red line, the white line marks the condition griQ — A, 
corresponding to the metastability of the upper minimum. Be- 
low the critical point, the free energy decreases when the dis- 
persion is tilted, indicating an instability. 



from A = 0. If the system is below the critical point, the 
free energy decreases, indicating an instability towards a 
Z2 symmetry broken state. The critical temperature in- 
creases when a non-zero tilt is chosen, as indicated by the 
red line. When the dispersion is tilted, the phase space 
density increases in the lower minimum, which results in a 
higher condensation temperature. For the non-interacting 
system the ratio of the critical temperatures for large tilt 
and no tilt is Tc,a^oo/^c,a = o = 2^^^. Furthermore, the 
condition | A| = griQ is shown by a white line. As discussed, 
this gives the order of magnitude of A for which the higher 
minimum of the free energy vanishes. This estimate is ac- 
curate for small temperatures, and gives an approximate 
energy scale for higher temperatures. We see that it is a 
large scale compared to the tilt energies studied in exper- 
iment. Therefore, the higher minimum is typically sta- 
ble, once the system is subcritical. In Fig. S12 we show 
the density fraction (tia — n^)/n^^ as function of the flux 
strength ^ for different temperatures. For high tempera- 
tures, a linear response to the tilt can be seen, while the 
appearance of a discontinuity for subcritical temperatures 
is indicative of a phase transition. 



ESTIMATION OF THE CRITICAL 
TEMPERATURE 

One can obtain a reasonable estimate for the number 
of Bose-condensed atoms at higher temperatures by ap- 
proximating the system as non-interacting and neglecting 
the trap in the x?/-plane. In that case, the total number 
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FIG S12. Density imbalance as a function of the flux. 

The density imbalance (n^ — n^^jn^^ between the two minima 
in the dispersion is calculated for three different temperatures 
as a function of the flux strength $. For subcritical temper- 
atures of the Z2 transition, the density imbalance has a dis- 
continuity at $ = TT. For supercritical temperatures we see a 
linear dependency. 



of atoms is given by the Bose statistics for the dispersion 
relation £:(k) together with the discrete level structure of 
the harmonic trap in z direction, 



A^Tot 



00 



£(fccc,fcy) + fe^z(r^z+l/2)- 



1 



) '. (S28) 



The Bose-condensate is, for given chemical potential /i, 
identified with the atom number A/q in the minimum of 
the dispersion relation. To work at fixed particle number, 
we adjust /i until the density equals a desired value, for 
which we use a realistic number (A^Tot/^Sites = 90 in a 
rhombic lattice of A^sites = 17 x 17). Further, we assume 
the experimental value hw^j J^^^^ = 3, and we use the ef- 
fective tunnelings for a given shaking modulation S (see 
Fig. SI). We normalise all quantities to the magnitude 
of the resulting tunneling matrix elements J®^. The re- 
sults for the number of excited atoms in the exact band 
structure are shown in Fig. S13a. The condensate is much 
stronger depleted at fluxes close to tt, indicating a stronger 
loss of /7(1) long-range order and a lower critical point for 
Bose condensation. This results in the same cusp-like be- 
haviour as obtained with the calculations in the weakly- 
interacting system with a harmonic approximation to the 
band structure (Fig. 5b of the main text). 

From the occupation at different k modes, one can also 
compute the peak width of the momentum distribution, 
similar to what is extracted from the experimental time- 
of- flight images and plotted in Fig. 5a of the main text. 
The result is shown in Fig. S13b. It reproduces well the 
qualitative behaviour of the experimental findings: the 
momentum peak gets broader closer to tt flux and at 
higher temperatures, pointing at a decrease of the U{1) 
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FIG S13. Calculations for the exact band structure in 
a non-interacting approximation. Data for A^Tot/A^sites — 

90 and A^sites = 17 x 17. a, The density of excited atoms 
riexc = Nexc/NTot is strongly enhanced at a flux strength of 
TT, resulting in a cusp, b, The FWHM of the momentum dis- 
tribution increases close to $ = tt and at higher temperatures, 
pointing at a decrease of U{1) long-range order. The values 
are normalised to the lowest result. 



long-range order. 



MONTE-CARLO SIMULATION 

In order to study the equilibrium states of the three 
dimensional ultracold ensemble, we simulate the system 
using classical Monte-Carlo (Metropolis algorithm). We 
start with a system of 42 x 75 tubes, where each tube 
is discretised into 33 sites using a discretisation length of 
az = I /im. This introduces an additional effective tunnel- 
ing term in z-direction J^ = h'^/{2ma1) = 10.7 J®^. Fur- 
thermore, the on-site repulsive interaction constant has to 
be rescaled to UsHe = 9iu/^z- 

The system is initialised using one of the two-fold de- 
generate classical ground states at zero flux with an initial 
total number of atoms of 2 x 10^. To ensure that the total 
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number of atoms remains constant, we set the chemical 
potential to /i = 230 J®^. We then perform single-site 
updates, where sites are chosen randomly and changes 
to the real and imaginary parts of the wavefunction are 
generated by sampling from a normal distribution whose 
width is adjusted to approach a step acceptance rate of 
roughly one half. After a thermalisation process that, de- 
pending on the temperature, consists in 10^ — 10^ Monte- 
Carlo steps per site (MCS), we start taking 100 snapshots 
of the systems. Between subsequent samples we perform 
sufficiently many MCS for both samples being completely 
uncorrelated: correlation is eradicated in under ten MCS, 
whereas our sampling frequency is between 100 and 1000 
MCS. For each sample we compute the chirality which is 
given by the relative visibility of the integrals over both 
interference peaks in quasimomentum space using a tri- 
angular mask. We repeat this process for flux values 
G [O.SItt, I.IQtt] and obtain the chirality as a function of 
the flux. The results are shown in Fig. 4. 



FIG S14. Monte-Carlo sample snapshot. Sample snap- 
shot of thermal equilibrium ensemble for flux strength $ = tt 
at T = 37 J^^ (so T > Tising), where the x^y-plane is cutting the 
tubes at the center of the trap: triangular plaquettes with neg- 
ative (positive) bosonic currents in red (blue). Higher (lower) 
color intensity represents higher (lower) absolute values of the 
bosonic current. In addition, regions with lower density are 
covered in white haze. Inset: Sample from center with arrows 
representing the phase at each lattice site. 
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